function [ft] = compare_stats(a, b, pflag, title, fid)
% compare_stats - provide comparison statistics between
% two groups of data - basic t test (unequal variances) and f test
% for ratio of variances. Optionally print the results to the screen
% fid is the file hande for output = 1 is stdout, 2 stderr, anything else is a file
%
% 12/9/98 P. Manis
%
[ft]= [];
if(isempty(a.x) | isempty(b.x))
   return;
end
% first use KS test to find if data is approximately normal in distribution,
% by comparing data with a normal distribution computed from the estimated
% mean and variance. Computation is taken from Press et al.
%
[ap] = ksone(a.x);
[bp] = ksone(b.x);
% check to see if it failed?
if(ap < 0.05 | bp < 0.05) % if either data set is not normally distributed
   [p] = wilcox_ranksum(a.x, b.x, 0.05); % run a wilcoxon rank-sums test.
   if(nargin > 2 & pflag)
      if(length(title) > 0)
         fprintf(fid, '%s: \n', title);
      end
		fprintf(fid, 'Data Not Normally distributed: A: p<%.3f  B: p<%.3f\n', ap, bp);
      fprintf(fid, 'Wilcoxon Rank Sums test:                          p = %.5f  %s\n', p, psig(p));
   end
else
   % compute standard parametric statistics using our tools (press et al versions)
   [f, fp, df1, df2] = f_test(a.x, b.x); % f test (probably should do first...)
   ft.f=f; ft.fp = fp; ft.df1 = df1; ft.df2 = df2;
   
   if(ft.fp < 0.05) % then we need to do the unequal variance version
      [t, tp, tdf] = tu_test(a.x, b.x); % unpaired heteroscedastic t test
      ft.t = t; ft.tp = tp; ft.tdf = tdf; ft.var=1;
   else
      [t, tp, tdf] = t_test(a.x, b.x); % unpaired equal variance t test
      ft.t = t; ft.tp = tp; ft.tdf = tdf; ft.var=0;
   end
   if(nargin > 2 & pflag)
      if(length(title) > 0)
         fprintf(fid, '%s:\n', title);
      end
         fprintf(fid, 'F-test: f = %.3f  df1 = %3d df2 = %3d             p = %.4f\n', ...
         ft.f, ft.df1, ft.df2, ft.fp);
      if(ft.var == 0)
         fprintf(fid, 'T-test (equal variance): t = %.3f  df = %3d      p = %.4f  %s\n', ...
            ft.t, ft.tdf, ft.tp, psig(ft.tp));
      else
         fprintf(fid, 'T-test (unequal variance): t = %.3f  df = %.1f  p = %.4f  %s\n', ...
            ft.t, ft.tdf, ft.tp, psig(ft.tp));
      end
      
   end
end

function [stars] = psig(p)
stars = ' ';
if(p <= 0.05 & p > 0.01)
stars = '*';
end
if(p <= 0.01 & p > 0.001)
   stars = '**';
end
if(p <= 0.001)
   stars = '***';
end




